rm(list = ls())
## llamando a las librerias que (creemos que) vamos a usar
options("install.lock"=FALSE)
if (!require ("pacman")) install.packages ("pacman")

pacman::p_load (dplyr
                , haven
                , ggplot2
                , ggpubr
                , ggthemes  ## si quieren mas themes
                , readstata13
                , readxl
                , sf
                , tidyverse
                , tidyr
                , units
                , viridis ## paleta de colores Viridis
                , wesanderson## p/usar paleta de colores de Wes Anderson
                , stringr
                , RColorBrewer
                , patchwork
                , Rmisc
                , lfe
                , stargazer
                , AER
                , haven
                , skimr
                , modelsummary
                , terra
                , fixest
                , vtable
                , did
                , cowplot
                , grid)

setwd("/Users/mateovasquez/Dropbox/")

base <- read.dta13("master_violence_landtitle_caracteristicas.dta")

#//3. Diff-in-Diff (Benchmark Specification) - Disaggregated 

names(base)

base <- base %>% 
  dplyr::mutate(total_violence = rowSums(across(c(sumattacks_public, sumattacks_guerrilla,sumattacks_paramilitary)), na.rm = TRUE),
                health_inst = rowSums(across(c(health_post, health_center,hospital)), na.rm = TRUE))


#Renaming Variables 

base <- base %>% 
  dplyr::rename("indian_title_area_cs" = "area_resguardos_indigenas_cs",
                "indian_title_area" = "area_resguardos_indigenas",
                "afro_title_area_cs" = "area_titul_comunidades_negras_cs",
                "afro_title_area" = "area_titul_comunidades_negras")

#Treatment 	

base <- base %>% 
  dplyr::mutate(post = ifelse(Año>1995,1,0),
                any_titleXpost = any_title*post,
                any_title2 = any_title,
                post2 = post)


#the mean level of violence in neighboring municipalities (1988–1995)
Neighbors <- base %>% 
  dplyr::filter(Año<=1995) %>% 
  dplyr::group_by(codmpio) %>% 
  dplyr::summarise(spatlagguerrsum100_100_1995 = mean(spatlagguerrsum100_100, na.rm = T))

base <- left_join(base,Neighbors,
                  by = "codmpio")

#Controls X post
base <- base %>%
  dplyr::mutate(
    altura_post = altura*post, #ELEVATION
    rainfall_post = rainfall*post, #RAINFALL
    coca_0_1_post = coca_0_1*post, # COCA
    presence_mines1560_post = presence_mines1560*post, #GOLD MINES
    num_encomiendas1560_post = num_encomiendas1560*post, #ENCOMIENDAS
    royalroaddist_post = royalroaddist*post, #distance to royal roads
    distancia_mercado_post = distancia_mercado*post, #nearest market towns (km)
    pop95_post = pop95*post, #municipal population in 1995 
    pct_minority85_post = pct_minority85*post, #share of minority population in 1985 (%)
    conflictos_1901_1917_post = conflictos_1901_1917*post,#historical conflict using data on the inci- dence of land conflict between 1901 and 1931
    conflictos_1918_1931_post = conflictos_1918_1931*post,#historical conflict using data on the inci- dence of land conflict between 1901 and 1931
    Violencia_48_a_53_post = Violencia_48_a_53*post, #presence of armed combat during La Violencia (1948–1953)
    anucraids1971_78_post = anucraids1971_78*post, #frequency of land invasions between 1971 and 1978
    spatlagguerrsum100_100_1995_post = spatlagguerrsum100_100_1995*post, #the mean level of violence in neighboring municipalities (1988–1995)
    ltotal_plots_rfmd_1960_85_post = ltotal_plots_rfmd_1960_85*post, # total number of plots reformed between 1960 and 1985 to proxy for rural grievances.
    totviolencia85_post = totviolencia85*post) #ojo: La Violencia (1948 - 1958)

# Applying asinh transformation to each specified variable directly in mutate and other variables

base <- base %>%
  mutate(
    as_total_violence = asinh(total_violence),
    as_sumattacks_public = asinh(sumattacks_public),
    as_sumattacks_paramilitary = asinh(sumattacks_paramilitary),
    as_sumattacks_guerrilla = asinh(sumattacks_guerrilla))

twfe.est_peas <- feols(as_sumattacks_paramilitary ~ i(Año,factor(any_title_placebo), 1995) | 
                    codmpio,  
                  data = base %>% dplyr::filter(gpacifica==1), cluster = "codmpio")

twfe.est_ind <- feols(as_sumattacks_paramilitary ~ i(Año,factor(any_title_reserve), 1995) | 
                    codmpio,  
                  data = base %>% dplyr::filter(gpacifica==1), cluster = "codmpio")

twfe.est.output_peas <- cbind.data.frame(twfe.est_peas$coefficients,
                                    SE = twfe.est_peas$se) %>% 
  rownames_to_column()

twfe.est.output_ind <- cbind.data.frame(twfe.est_ind$coefficients,
                                         SE = twfe.est_ind$se) %>% 
  rownames_to_column()

twfe.est.output_peas <- as.data.frame(twfe.est.output_peas)
twfe.est.output_ind <- as.data.frame(twfe.est.output_ind)

twfe.est.output_peas$Year <- substr(twfe.est.output_peas$rowname, 6, 9)
twfe.est.output_peas$Group <- substr(twfe.est.output_peas$rowname, 38, 38)

twfe.est.output_ind$Year <- substr(twfe.est.output_ind$rowname, 6, 9)
twfe.est.output_ind$Group <- substr(twfe.est.output_ind$rowname, 38, 38)

twfe.est.output_peas <- twfe.est.output_peas %>% 
  dplyr::rename("Estimate" = "twfe.est_peas$coefficients") %>% 
  dplyr::filter(Group == 1) %>% 
  dplyr::mutate(Year = as.numeric(Year),
                LB_CI = Estimate-(1.96*SE),
                UB_CI = Estimate+(1.96*SE),
                Model = "TWFE",
                interval = ifelse(Year < 1985, 1980,
                                  ifelse(Year < 1989,1985,
                                         ifelse(Year < 1994,1990,
                                                ifelse(Year < 2000,1995,
                                                       ifelse(Year < 2005,2000,
                                                              ifelse(Year < 2010,2005,2010))))))) %>% 
  dplyr::filter(!is.na(Year)) %>% 
  dplyr::select(-c(SE,rowname))

twfe.est.output_ind <- twfe.est.output_ind %>% 
  dplyr::rename("Estimate" = "twfe.est_ind$coefficients") %>% 
  dplyr::filter(Group == 1) %>% 
  dplyr::mutate(Year = as.numeric(Year),
                LB_CI = Estimate-(1.96*SE),
                UB_CI = Estimate+(1.96*SE),
                Model = "TWFE",
                interval = ifelse(Year < 1985, 1980,
                                  ifelse(Year < 1989,1985,
                                         ifelse(Year < 1994,1990,
                                                ifelse(Year < 2000,1995,
                                                       ifelse(Year < 2005,2000,
                                                              ifelse(Year < 2010,2005,2010))))))) %>% 
  dplyr::filter(!is.na(Year)) %>% 
  dplyr::select(-c(SE,rowname))

results_peas <- twfe.est.output_peas %>%
  dplyr::group_by(interval, Group) %>%
  dplyr::summarise(Estimate_mean = mean(Estimate, na.rm = TRUE),
                   LB_CI_mean = mean(LB_CI, na.rm = TRUE),
                   UB_CI_mean = mean(UB_CI, na.rm = TRUE))

results_ind <- twfe.est.output_ind %>%
  dplyr::group_by(interval, Group) %>%
  dplyr::summarise(Estimate_mean = mean(Estimate, na.rm = TRUE),
                   LB_CI_mean = mean(LB_CI, na.rm = TRUE),
                   UB_CI_mean = mean(UB_CI, na.rm = TRUE))

# Assuming 'results' is your data frame and 'Group' is the column you're adjusting
results_peas$Group <- "Campesino"
results_ind$Group <- "Indigenous"

results <- rbind(results_peas,
                 results_ind)
####

twfe.est <- feols(as_sumattacks_paramilitary ~ i(Año,factor(any_title2), 1995) | 
                    codmpio,  
                  data = base %>% dplyr::filter(gpacifica == 1), cluster = "codmpio")

twfe.est.output <- cbind.data.frame(twfe.est$coefficients,
                                    SE = twfe.est$se) %>% 
  rownames_to_column()

twfe.est.output <- as.data.frame(twfe.est.output)

twfe.est.output$Year <- substr(twfe.est.output$rowname, 6, 9)
twfe.est.output$Group <- substr(twfe.est.output$rowname, 31, 31)

twfe.est.output <- twfe.est.output %>% 
  dplyr::rename("Estimate" = "twfe.est$coefficients") %>% 
  dplyr::mutate(Year = as.numeric(Year),
                LB_CI = Estimate-(1.96*SE),
                UB_CI = Estimate+(1.96*SE),
                Model = "TWFE",
                interval = ifelse(Year < 1985, 1980,
                                  ifelse(Year < 1989,1985,
                                         ifelse(Year < 1994,1990,
                                                ifelse(Year < 2000,1995,
                                                       ifelse(Year < 2005,2000,
                                                              ifelse(Year < 2010,2005,2010))))))) %>% 
  dplyr::filter(!is.na(Year)) %>% 
  dplyr::select(-c(SE,rowname))


results_titled <- twfe.est.output %>%
  dplyr::group_by(interval, Group) %>%
  dplyr::summarise(Estimate_mean = mean(Estimate, na.rm = TRUE),
                   LB_CI_mean = mean(LB_CI, na.rm = TRUE),
                   UB_CI_mean = mean(UB_CI, na.rm = TRUE))

# Assuming 'results' is your data frame and 'Group' is the column you're adjusting

results_titled <- results_titled %>% 
  dplyr::filter(Group == 1)

results_titled$Group <- "Titled"


results <- rbind(results_peas,
                 results_ind,
                 results_titled)

# Then, your plotting code follows

Paramilitary_pac_camp<- ggplot(results %>% dplyr::filter(Group == "Titled" |
                                                                   Group == "Campesino"), aes(x = interval, y = Estimate_mean, shape = as.factor(Group), color = as.factor(Group))) +
  geom_line(aes(linetype=as.factor(Group)), size = 1) +
  scale_linetype_manual(values=c("solid", "dashed"),
                        breaks=c("Titled", "Campesino"),
                        labels=c("Titled (Law 70)", "Campesino")) +
  geom_point(size = 3) +
  scale_shape_manual(values = c("Titled" = 18, "Campesino" = 15),
                     breaks=c("Titled", "Campesino"),
                     labels=c("Titled (Law 70)", "Campesino")) +
  scale_color_manual(values = c("Titled" = "#000033", "Campesino" = "#990000"),
                     breaks=c("Titled", "Campesino"),
                     labels=c("Titled (Law 70)", "Campesino")) + # Define colors for the groups
  geom_errorbar(aes(ymin=LB_CI_mean, ymax=UB_CI_mean, group= factor(Group)), size=0.2, width=0.2) +
  labs(x = "Interval", y = "Marginal Effect Coefficient", title = "Titled (Law 70) Vs Campesino") +
  geom_hline(yintercept = mean(base$as_sumattacks_paramilitary[base$gpacifica==1]), col = "red") +
  geom_vline(xintercept = 1995, col = "black") +
  theme_calc() +
  theme(legend.position = c(0.1, 0.9),
        legend.justification = c(0, 1),
        legend.text = element_text(size = 12),
        legend.key.size = unit(1.5, "lines"),
        legend.background = element_rect(fill = "white", colour = "black")) +
  scale_x_continuous(breaks = seq(1980, 2010, by = 5), limits = c(1980, 2010)) +
  scale_y_continuous(breaks = seq(-.5, 1, by = .5), limits = c(-.5, 1)) +
  guides(shape = guide_legend(title = NULL), linetype = guide_legend(title = NULL), color = guide_legend(title = NULL))

# Then, your plotting code follows

Paramilitary_pac_Peasantcamp <- ggplot(results %>% dplyr::filter(Group == "Campesino" |
                                                                   Group == "Indigenous"), aes(x = interval, y = Estimate_mean, shape = as.factor(Group), color = as.factor(Group))) +
  geom_line(aes(linetype=as.factor(Group)), size = 1) +
  scale_linetype_manual(values=c("solid", "dashed"),
                        breaks=c("Campesino", "Indigenous"),
                        labels=c("Campesino", "Indigenous")) +
  geom_point(size = 3) +
  scale_shape_manual(values = c("Campesino" = 18, "Indigenous" = 15),
                     breaks=c("Campesino", "Indigenous"),
                     labels=c("Campesino", "Indigenous")) +
  scale_color_manual(values = c("Campesino" = "#336600", "Indigenous" = "orange"),
                     breaks=c("Campesino", "Indigenous"),
                     labels=c("Campesino", "Indigenous")) + # Define colors for the groups
  geom_errorbar(aes(ymin=LB_CI_mean, ymax=UB_CI_mean, group= factor(Group)), size=0.2, width=0.2) +
  labs(x = "Interval", y = "Marginal Effect Coefficient", title = "Campesino Vs Indigenous") +
  geom_hline(yintercept = mean(base$as_sumattacks_paramilitary[base$gpacifica==1]), col = "red") +
  geom_vline(xintercept = 1995, col = "black") +
  theme_calc() +
  theme(legend.position = c(0.1, 0.9),
        legend.justification = c(0, 1),
        legend.text = element_text(size = 12),
        legend.key.size = unit(1.5, "lines"),
        legend.background = element_rect(fill = "white", colour = "black")) +
  scale_x_continuous(breaks = seq(1980, 2010, by = 5), limits = c(1980, 2010)) +
  scale_y_continuous(breaks = seq(-.5, 1, by = .5), limits = c(-.5, 1)) +
  guides(shape = guide_legend(title = NULL), linetype = guide_legend(title = NULL), color = guide_legend(title = NULL))

####

Paramilitary_pac_ind <- ggplot(results %>% dplyr::filter(Group == "Titled" |
                                                      Group == "Indigenous"), aes(x = interval, y = Estimate_mean, shape = as.factor(Group), color = as.factor(Group))) +
  geom_line(aes(linetype=as.factor(Group)), size = 1) +
  scale_linetype_manual(values=c("solid", "dashed"),
                        breaks=c("Titled", "Indigenous"),
                        labels=c("Titled (Law 70)", "Indigenous")) +
  geom_point(size = 3) +
  scale_shape_manual(values = c("Titled" = 18, "Indigenous" = 15),
                     breaks=c("Titled", "Indigenous"),
                     labels=c("Titled (Law 70)", "Indigenous")) +
  scale_color_manual(values = c("Titled" = "#000033", "Indigenous" = "#990000"),
                     breaks=c("Titled", "Indigenous"),
                     labels=c("Titled (Law 70)", "Indigenous")) + # Define colors for the groups
  geom_errorbar(aes(ymin=LB_CI_mean, ymax=UB_CI_mean, group= factor(Group)), size=0.2, width=0.2) +
  labs(x = "Interval", y = "Marginal Effect Coefficient", title = "Law 70 Vs Indigenous") +
  geom_hline(yintercept = mean(base$as_sumattacks_paramilitary[base$gpacifica==1]), col = "red") +
  geom_vline(xintercept = 1995, col = "black") +
  theme_calc() +
  theme(legend.position = c(0.1, 0.9),
        legend.justification = c(0, 1),
        legend.text = element_text(size = 12),
        legend.key.size = unit(1.5, "lines"),
        legend.background = element_rect(fill = "white", colour = "black")) +
  scale_x_continuous(breaks = seq(1980, 2010, by = 5), limits = c(1980, 2010)) +
  scale_y_continuous(breaks = seq(-.5, 1, by = .5), limits = c(-.5, 1)) +
  guides(shape = guide_legend(title = NULL), linetype = guide_legend(title = NULL), color = guide_legend(title = NULL))


Main_plot <- ggarrange(Paramilitary_pac_camp, Paramilitary_pac_ind, Paramilitary_pac_Peasantcamp, ncol = 3, nrow = 1)

annotate_figure(Main_plot, top = text_grob("Paramilitary Attacks", 
                                           color = "black", face = "bold", size = 14))

